﻿#define _CRT_SECURE_NO_WARNINGS
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include "Interp.h"
#include "RBtree.h"

inline double Interp::GetMin(double val1, double val2)
{
	return ((val1) < (val2) ? (val1) : (val2));
}

inline double Interp::GetMax(double val1, double val2)
{
	return ((val1) > (val2) ? (val1) : (val2));
}

inline int Interp::POWMINUS(int x)
{
	return (1 - 2 * (x % 2));
}

inline double* Interp::POINTELEM(double* BASE, unsigned int N, unsigned int X, unsigned int Y)
{
	return (BASE + (X - 1) * N + Y - 1);
}

bool Interp::MatLoadCSV(const char* path, double*& data, unsigned int& rows, unsigned int& cols)
{
	FILE* file;
	file = fopen(path, "r");
	if (file == NULL)
	{
		printf("读取CSV文件时发生错误：无法打开文件\n");
		return false;
	}
	rows = 0;
	int c0 = -1, c1;
	//char buffer[1024*5];
	char* buffer = (char*)malloc(MAXROW * sizeof(char));
	while (fgets(buffer, MAXROW * sizeof(char), file) != NULL)
	{
		rows++;
		c1 = 0;
		for (int i = 0; i < MAXROW * sizeof(char); i++)
		{
			if (buffer[i] == '\r' || buffer[i] == '\n' || buffer[i] == '\0')
				break;
			if (buffer[i] == ',')
				c1++;
		}
		if (c0 != -1 && c0 != c1)
		{
			printf("读取CSV文件时发生错误：每行数据列数不一致\n");
			return false;
		}
		c0 = c1;
	}
	free(buffer);
	cols = c0 + 1;
	data = (double*)malloc(rows * cols * sizeof(double));
	rewind(file);
	for (int i = 1; i <= rows; i++)
		for (int j = 1; j <= cols; j++)
			if (j == cols)
				fscanf(file, "%lf", POINTELEM(data, cols, i, j));
			else
				fscanf(file, "%lf,", POINTELEM(data, cols, i, j));
	fclose(file);
	return true;
}

bool Interp::MatLoadTXT(const char* path, double*& data, unsigned int& rows, unsigned int& cols)
{
	FILE* file;
	file = fopen(path, "r");
	if (file == NULL)
	{
		printf("读取TXT文件时发生错误：无法打开文件\n");
		return false;
	}
	rows = 0;
	int c0 = -1, c1;
	char* buffer = (char*)malloc(MAXROW * sizeof(char));
	while (fgets(buffer, MAXROW * sizeof(char), file) != NULL)
	{
		rows++;
		c1 = 0;
		for (int i = 0; i < MAXROW * sizeof(char) - 1; i++)
		{
			if (buffer[i + 1] == ' ' || buffer[i + 1] == '\t')
			{
				if (buffer[i] != ' ' && buffer[i] != '\t')
					c1++;
			}
			else if (buffer[i + 1] == '\r' || buffer[i + 1] == '\n' || buffer[i + 1] == '\0')
			{
				if (buffer[i] != ' ' && buffer[i] != '\t')
					c1++;
				break;
			}
		}
		free(buffer);
		if (c0 != -1 && c0 != c1)
		{
			printf("读取TXT文件时发生错误：每行数据列数不一致\n");
			return false;
		}
		c0 = c1;
	}
	cols = c0;
	data = (double*)malloc(rows * cols * sizeof(double));
	rewind(file);
	for (int i = 1; i <= rows; i++)
		for (int j = 1; j <= cols; j++)
			fscanf(file, "%lf", POINTELEM(data, cols, i, j));
	fclose(file);
	return true;
}

Interp::Interp()
{
	OnInti = false;
}
void Interp::Delete()
{
	if (OnInti) 
	{
		OnInti = false;
		for (int i = 0;i < dim;i++)
			free(argu[i]);
		free(data);
	}
}
void Interp::Print()
{
	printf("%d\n", dim);
	for (int i = 0;i < dim;i++)
		printf("%d ", num[i]);
	printf("\n");
	for (int i = 0;i < dim;i++)
	{
		for (int j = 0;j < num[i];j++)
			printf("%lf ", argu[i][j]);
		printf("\n");
	}
	for (int i = 0;i < tol;i++)
		printf("%lf ", data[i]);
}
bool Interp::LoadBIN(const char* path) 
{
	if (OnInti)
		Delete();
	OnInti = true;
	FILE* file;
	file = fopen(path, "rb");
	if (file == NULL)
	{
		printf("读取二进制文件时发生错误：无法打开文件\n");
		OnInti = false;
		return false;
	}
	fread(&dim, sizeof(unsigned int), 1, file);
	tol = 1;
	for (int i = 0;i < dim;i++) 
	{
		fread(&num[i], sizeof(unsigned int), 1, file);
		tol *= num[i];
	}
	for (int i = 0;i < dim;i++) 
	{
		argu[i] = (double*)malloc(num[i] * sizeof(double));
		for (int j = 0;j < num[i];j++)
			fread(&argu[i][j], sizeof(double), 1, file);
	}
	for (int i = 0; i < dim; i++)
		OnEq[i] = true;
	for (int i = 0;i < dim;i++) 
		for (int j = 1;j < num[i] - 1;j++)
			if (argu[i][j - 1] + argu[i][j + 1] != 2 * argu[i][j])
			{
				OnEq[i] = false;
				break;
			}
	data = (double*)malloc(tol * sizeof(double));
	for (int i = 0;i < tol;i++)
		fread(&data[i], sizeof(double), 1, file);
	fclose(file);
	return true;
}
bool Interp::SaveBIN(const char* path) 
{
	if (!OnInti)
	{
		printf("保存二进制文件时发生错误：插值表未初始化\n");
		return false;
	}
	FILE* file;
	file = fopen(path, "wb");
	if (file == NULL)
	{
		printf("保存二进制文件时发生错误：无法打开文件\n");
		return false;
	}
	fwrite(&dim, sizeof(unsigned int), 1, file);
	for (int i = 0;i < dim;i++)
		fwrite(&num[i], sizeof(unsigned int), 1, file);
	for (int i = 0;i < dim;i++)
		for (int j = 0;j < num[i];j++)
			fwrite(&argu[i][j], sizeof(double), 1, file);
	for (int i = 0;i < tol;i++)
		fwrite(&data[i], sizeof(double), 1, file);
	fclose(file);
	return true;
}
double* Interp::PointData(unsigned int* index) 
{
	if (!OnInti)
	{
		printf("求指定插值点数据的地址时发生错误：插值表未初始化\n");
		return NULL;
	}
	unsigned int sum = index[0];
	for (int i = 1;i < dim;i++)
		sum = sum * num[i] + index[i];
	return data + sum;
}
void Interp::GetIndex(const double* x, unsigned int* index, bool term)
{
	double val = 0;
	for (int i = 0;i < dim;i++) 
	{
		val = GetMax(x[i], argu[i][0]);
		val = GetMin(val, argu[i][num[i] - 1] - 1e-5 * term);
		if (OnEq[i])
			index[i] = floor((val - argu[i][0]) / (argu[i][1] - argu[i][0]));
		else 
		{
			unsigned int low = 0;
			unsigned int high = num[i] - 1;
			unsigned int middle = 0;
			while (low < high) 
			{
				middle = (low + high) / 2;
				if (val == argu[i][middle])
					break;
				if (val < argu[i][middle])
					high = middle - 1;
				else
					low = middle + 1;
			}
			if (val == argu[i][middle])
				index[i] = middle;
			else if (val < argu[i][low])
				index[i] = low - 1;
			else
				index[i] = low;
		}
	}
}
double Interp::GetVal(const double* x, int t) 
{
	if (t == dim) //递归完成，开始回溯
	{
		unsigned int index0[MAXNUM];
		unsigned int index[MAXNUM];
		GetIndex(x,index0);
		double nume = 1;
		double deno = 1;
		double val = 0;
		for (int i = 0;i < dim;i++) 
		{
			val = GetMax(x[i], argu[i][0]);
			val = GetMin(val, argu[i][num[i] - 1] - 1e-5);
			index[i] = index0[i] + recu[i];
			nume *= POWMINUS(recu[i]) * (argu[i][index0[i] + 1 - recu[i]] - val);
			deno *= argu[i][index0[i] + 1] - argu[i][index0[i]];
		}
		return nume * (*PointData(index)) / deno;
	}
	double sum = 0;
	recu[t] = 0; //左分支
	sum += GetVal(x, t + 1);
	recu[t] = 1; //右分支
	sum += GetVal(x, t + 1);
	return sum;
}
bool Interp::Load(const char* path, int _dim, bool fillZeros)
{
	if (OnInti)
		Delete();
	OnInti = true;
	dim = _dim;
	double* mdata;
	unsigned int rows, cols;
	char ext[4];
	strncpy(ext, path + strlen(path) - 3, 4);
	if (!strcmp(ext, "csv") || !strcmp(ext, "CSV")) 
	{
		MatLoadCSV(path, mdata, rows, cols);
	}
	else if (!strcmp(ext, "txt") || !strcmp(ext, "TXT")) 
	{
		MatLoadTXT(path, mdata, rows, cols);
	}
	else 
	{
		printf("处理插值表时发生错误：只支持.csv文件和.txt文件\n");
		return false;
	}
	rbtree* T = rbtree_create();
	double* temp = (double*)malloc(MAXPOINTNUM * sizeof(double));
	tol = 1;
	for (int i = 0;i < dim;i++) 
	{
		for (int j = 1; j <= rows; j++)
			rbtree_insert(T, *POINTELEM(mdata, cols, j, i + 1));
		unsigned int index = 0;
		rbtree_out(T, T->root, temp, index);
		num[i] = index;
		tol *= num[i];
		argu[i] = (double*)malloc(num[i] * sizeof(double));
		memcpy(argu[i], temp, num[i] * sizeof(double));
		rbtree_destroy(T, T->root);
		T->root = T->nil;
	}
	free(T->nil);
	free(temp);
	for (int i = 0; i < dim; i++)
		OnEq[i] = true;
	for (int i = 0;i < dim;i++)
		for (int j = 1;j < num[i] - 1;j++)
			if (argu[i][j - 1] + argu[i][j + 1] != 2 * argu[i][j])
			{
				OnEq[i] = false;
				break;
			}
	if (rows < tol && !fillZeros)
	{
		printf("处理插值表时发生错误：插值点个数过少\n");
		free(mdata);
		for (int i = 0;i < dim;i++)
			free(argu[i]);
		OnInti = false;
		return false;
	}
	else if (rows > tol && !fillZeros)
	{
		printf("处理插值表时发生错误：插值点个数过多\n");
		free(mdata);
		for (int i = 0;i < dim;i++)
			free(argu[i]);
		OnInti = false;
		return false;
	}
	data = (double*)malloc(tol * sizeof(double));
	if (fillZeros)
		memset(data, 0, tol * sizeof(double));
	int* flag = (int*)malloc(tol * sizeof(int));
	memset(flag, 0, tol * sizeof(int));
	for (int i = 1;i <= rows;i++) 
	{
		double x[MAXNUM];
		for (int j = 0; j < dim; j++)
			x[j] = *POINTELEM(mdata, cols, i, j + 1);
		unsigned int index[MAXNUM];
		GetIndex(x, index, false);
		double* pos = PointData(index);
		*pos = *POINTELEM(mdata, cols, i, dim + 1);
		int* pflag = flag + (unsigned int)(pos - data);
		if (*pflag) 
		{
			printf("处理插值表时发生错误：%d行数据与%d行数据发生冲突\n", i, *pflag);
			free(flag);
			Delete();
			free(mdata);
			OnInti = false;
			return false;
		}
		*pflag = i;
	}
	bool check = true;
	for (int i = 0;i < tol;i++) 
		if (!flag[i]) 
		{
			check = false;
			break;
		}

	if (!check && !fillZeros) 
	{
		printf("处理插值表时发生错误：插值点缺失\n");
		free(flag);
		Delete();
		free(mdata);
		OnInti = false;
		return false;
	}
	free(flag);
	free(mdata);
	return true;
}
